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Abstract 

The paper presents new results on localisation and transmission of flexural waves 
in a structured plate containing a semi-infinite two-dimensional array of rigid pins. In 
particular, surface waves are identified and studied at the interface boundary between 
the homogeneous part of the flexural plate and the part occupied by rigid pins. A 
formal connection has been made with the dispersion properties of flexural Bloch waves 
in an infinite doubly periodic array of rigid pins. Special attention is given to regimes 
corresponding to standing waves of different types as well as Dirac-like points, that may 
occur on the dispersion surfaces. A single half-grating problem, hitherto unreported 
in the literature, is also shown to bring interesting solutions. 

1 Introduction 

The advent of designer materials such as metamaterials, photonic crystals and micro-structured 
media that are able to generate effects unobtainable by natural media, such as Pendry’s flat 
lens m, is driving a revolution in materials science. Many of these ideas originate in elec¬ 
tromagnetism and optics, but are now percolating into other wave systems such as those of 
elasticity, acoustics or the idealised Kirchhoff-Love plate equations for flexural waves, with 
this analogue of photonic crystals being labelled as platonics [TH]. Many of the effects from 
photonics also appear in the flexural wave context albeit with some changes due to the 
biharmonic nature of the Kirchhoff-Love equation: ultra-refraction and negative refraction 
PIZ!, Dirac-like cones, Dirac-cone cloaking and related effects [alETlIIS] amongst others. 

The Kirchhoff-Love equations are good approximations within their realm of applica¬ 
bility (for instance with /i, A as plate thickness and typical wavelength respectively then 
/i/A <C 1 is required) and capture much of the essence of the wave physics, hence their 
emergent popularity. The system is relatively simple so analytic results for infinite periodic 
structured plates pinned, say, at regular points readily emerge [I6], with much of this earlier 
work reviewed in m; more recently multipole methods IZB, extending pins to cylinders, 
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or high-frequency homogenisation approaches [2] to get effective continuum equations that 
encapsulate the microstructure, have emerged. For finite pinned regions of a plate, a Green’s 
function approach [S] leads to rapid numerical solutions, or for an infinite grating, one may 
employ an elegant methodology for exploring Rayleigh-Bloch modes. This includes exten¬ 
sions to stacks of gratings and the trapping and filtering of waves [221 ESI [lU which further 
exemplify this approach. Problems such as these start to pose questions about semi-infinite 
gratings, or edge states in semi-infinite lattices, and our aim is to generate the relevant exact 
solutions. 

For semi-infinite cracks, and other situations where the change in boundary condition 
along a line is of interest, the classical Wiener-Hopf technique [SH ] for solving integral 
equations is highly developed and often used for mixed boundary value problems in contin¬ 
uum mechanics; indeed the seminal Wiener-Hopf treatise by [23] covers this aspect almost 
exclusively. However, lesser known and presented only as an exercise (4.10, p. 173-4 [H]) in 
[23], is its application to continuum discrete problems such as gratings. As such, it has 
seen notable application to the fracture of discrete lattice systems as reviewed by [23], to 
antenna design [30] and various semi-infinite grating and lattice scattering problems in acous¬ 
tics [ini 129] where the early work of [12] (hereafter referred to as HK), motivated by [13] , 
provides a wealth of useful information. 

The corresponding exact solutions for the platonic system are not available and, given 
the current interest in platonics and the versatility of these in, say, asymptotic schemes 
and the key insight given into the physics and results, we aim to provide these here. Of 
particular interest are regimes for which one observes the interplay between grating modes 
and Floquet-Bloch waves in a full doubly periodic infinite structure. In particular, these 
include frequencies and wave vector components corresponding to stationary points on the 
dispersion surfaces as well as neighbourhoods of the Dirac-like points. Additionally, we 
address a fundamental question as to whether a simple, yet surprisingly physically rich, 
structure such as a half-plane of rigid pins in a plate supports flexural interfacial modes; 
by which we mean waves that propagate along the interface of the platonic crystal and the 
homogeneous part of the biharmonic plate. 

Since the underlying mathematical structure of a semi-infinite linear grating and that of 
a semi-infinite lattice are very similar we choose to treat both examples within this paper. 
For clarity, we shall refer to a grating when we have a single half-line of scatterers (see figure 
[^a)) and to a lattice when there is a lattice of scatterers in a half-plane (see figure [^b)); 
throughout this paper we shall assume that we have point clamped scatterers although the 
analysis is readily extended to more complicated conditions holding at a point. In section 
we formulate the problem drawing upon the discrete Wiener-Hopf technique used by HK 
for the semi-infinite diffraction grating governed by a Helmholtz operator; much is directly 
applicable provided we replace the point source Green’s function for the two-dimensional 
Helmholtz operator with that for the biharmonic operator. In some regards, difficulties 
associated with the Helmholtz Green’s function’s logarithmic singularity are bypassed for 
the biharmonic case, where there is no singularity. We turn our attention to the linear 
grating first, in section m and follow that with an analysis of the lattice in section |2.2[ as 
noted above this is a parametrically and physically rich problem displaying a surprisingly 
wide range of behaviours and we outline methodologies for both interpreting and designing 
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Figure 1: (a) A single semi-infinite grating of rigid pins in a biharmonic plate with spacing 
s. The pins are represented by the large black discs. The forces/intensities are represented 
by the coefficients A*, and displacements, used in the Wiener-Hopf approach, by 6j. (b) The 
semi-infinite lattice of pins, consisting of an array of infinite gratings with period dy in the 
vertical direction and horizontal spacing d^- 

specific behaviours; we consider these in sections |3.1113.2| for the semi-infinite grating and 
lattice respectively. Concluding remarks are drawn together in section 


2 Formulation 

The methodology primarily adopted here, at least for the exact solution, is that of discrete 
Wiener-Hopf drawing upon the treatment of the Helmholtz semi-infinite grating by HK, 
which we briefly review here. To avoid issues with singularities in their problem, HK consider 
cylinders of finite radius, small compared to wavelength, and widely spaced. The diffracted 
field consists of a cylindrical wave plus a set of plane waves that possess amplitudes and 
directions of propagation identical to those for an infinite grating, but for the semi-infinite 
case, they do not exist everywhere. Instead the end-effects lead to shadow boundaries defined 
by lines drawn from the end of the grating to infinity along propagation directions of various 
plane waves. The cylindrical wave has sharp zeros followed by a maximum for certain 
directions; these zeros occur in the directions that would be shadow boundaries if the incident 
wave travels directly into the end of the grating - a class of resonant cases characterised by 
a diffracted wave that travels parallel to the grating. HK note that for the wave that travels 
parallel and out of the grating, the cylindrical wave and all the other plane waves vanish. 
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( 1 ) 


Assuming isotropic scatterers HK, see also na, pose the solution in the form 

oo 

U{r,e) = ^A^Hl,^\f5rr,), 

n=0 

for some coefficients to be determined and Hq^\z) being the usual Hankel function of 
the first kind with wavenumber /5, and representing the distance between the observation 
point (r, 9) and the nth scattering element on the a:-axis; for an infinite array these just diher 
by a phase factor and the analysis is simpler. The assumption that the grating’s elements are 
small, and therefore scatter isotropically, are not formally required for the pinned biharmonic 
plate as this is exactly true. 

In HK, the spacing between the scatterers is assumed to be large compared with the 
wavelength, and it must not be an integral or half-integral multiple of the wavelength. 
These assumptions are required by the mathematical analysis employed. The reason for the 
limitation on the integral or half-integral multiples of wavelength is to ensure that branch 
cuts for 2 : = exp (iPs) and 2 = exp (—iPs) are distinct (see appendix B for a discussion of 
branch cuts for the biharmonic plate) and s is the spacing of the pins. 

To determine the scattered wave coefficients An in 0 . the boundary conditions generate 
an infinite system of linear equations solvable using discrete Wiener-Hopf. In the neighbour¬ 
hood of the nth scatterer, the field consists of the incident plane wave, the sum of waves 
scattered to the nth element by the other elements and the wave scattered by the nth element 
itself: 

OO 

gin./3cosv- + An,H^^\p\n - m\s) + H^^\pa) = 0, n = 0,1, 2,.... (2) 

m=0 

m^n 

Here a is the radius of the individual elements. The term H^\pa) is treated separately 
because of the logarithmic singularity arising for the Hankel function at the origin. Although 
a small radius a is assumed by HK, the singularity poses problems as a —?• 0 and much of 
HK involves dealing with this; for the biharmonic case we discuss in this article, the Green’s 
function is finite at the source point rather than diverging logarithmically. 

2.1 Semi-infinite grating of rigid pins in a biharmonic plate 

For the Kirchhoff-Love plate equation 

A^G{r)-P^G{v) = 5{r-r'), (3) 

with an array of pins, the solution, cf. 0 , is 

OO 

U{r,0) = Y^AnG{Ppn) (4) 

n=0 

with the An to be determined, and where the free-space Green’s function 

G(Ppn) = ^ (GHfiPn) + |a'o(/?P«)) (5) 
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is a function of position = |r — r^|, with = (ns, 0) for spacing s and n G Z. As in 0 
the Green’s function contains but this is now augmented by Kq which is the modified 
Bessel function [T]; this Green’s function is bounded at r = r'. 

We take rigid pins located at points (ns, 0) for n > 0 and to set up the system of equations 
we introduce displacements bn at (ns, 0 ) for all n; imposing zero displacements at the pins, 
= 0 for n > 0 but are unknown for n < 0. For the intensities An in @ since there are no 
sources on the left-hand side, = 0 Vn > 0. The analogue of a is 

oo 

gm*/3cosV- ^ AmG(^\n - m\s) =bn, 71 E Z. (6) 

m=0 

We now employ the z-transform, multiply ^ by z'^ for z complex, and sum over all n\ 

oo oo oo oo 

Y + Y ^'"^mG{(3\7l - 7n\s) = Z^bn- (7) 

n=—oo n=—oo m=0 n=—oo 

To transform ([^ into a single functional equation of the Wiener-Hopf type, it is convenient 
to define A{z) and a kernel function IC{z) as 


A(z) = 

171=0 

OO ^. 

j=-oo 


( 8 ) 

(9) 


Notably, for IC(z), the j = 0 term is constant, i/SP‘^, which replaces the awkward H^\l3a) 
term of ([^. Using these functions, 0 becomes 


F{z) 7 A{z)K:{z) = B{z), 


and the forcing function, in this case, is 


OO 

where B{z) = 

n=l 


F(z) 




is^ cos'ip\'^ 


( 10 ) 


( 11 ) 


but could take different forms if the forcing were altered. 

Equation (10) is the starting point for the Wiener-Hopf technique with key steps being 
the product factorization of the kernel function }C{z) into two factors which are analytic in 
given but different regions, and factorization of the forcing function. Before proceeding we 
turn to the semi-infinite lattice as the formulation is almost identical. 


2.2 Semi-infinite lattice of rigid pins in a biharmonic plate 

We now replace each pin in the semi-infinite grating with an infinite grating in the vertical 
direction, as illustrated in figure [^b). The properties of the infinite grating are well known 
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e.g. [ID] and the quasi-periodic grating Green’s function 


Gl{l3,X]Ky,dy) = 


j=-oo 


H, 


( 1 ) 


P\IUds)^ +I + ^A’o ( l}\J{idyf + x‘ 


( 12 ) 

plays an important role. Here dy is the vertical period (in section 2.1 we used s to denote 
the period for the single grating) and Ky is the corresponding Bloch parameter. The gratings 
are separated by integer multiples of spacing d^ and the semi-infinite lattice is created from 
columns of infinite gratings (or rows of semi-infinite gratings). The analogue of ([^ emerges 
as 


e”''-’’”* + Xi AnGMn - m|4; Kj, d,) = 6„, n e Z. 


(13) 


m=0 


We repeat the procedure for the semi-infinite grating with the only difference being that the 
kernel function is now connected with the doubly quasi-periodic Green’s function: 


oo oo 


A{z))C{z) = yy AmGl_^{(3\n - m|4; Ky, d.. 




(14) 


n=—oo m=0 


when z = e 


it^x dx 


with Kx = P cos'll). 



Figure 2: The unit circle, 12+, 12_ and contours (7+ and G- which define the boundary of 
the ring of analyticity. The integrals on (7+ and (7_ are evaluated counter-clockwise and 
clockwise respectively. 


2.3 Wiener-Hopf 


The crux of the Wiener-Hopf methodology is to unravel (10) using the analyticity properties 
of the unknown functions A{z),B{z) and the known functions JC{z), F{z). To proceed we 
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define domains 0+ and Q- snch that 


t7+ = {z : \z\ < c+}, i7_ = {z : \z\ > c_}, 


( 15 ) 


where (7+ is a circle of radins c+ = 1 + 5 and C- is a circle of radius c_ = 1 — 5 for 
0 < 5 <C 1 (see figure ; the intersection 17+ fl 17_ describes an annulus of analyticity in 
the neighbourhood of the unit circle. From the definitions of A(z) Q and B{z) ( [lol ), and 
assuming appropriate decay at infinity, we identify them as + and — functions analytic in 
17+ and 17_ respectively. We then attempt to separate (10) into sides entirely analytic in 
17+; the only way both can be equal, given an overlapping region of analyticity, is for them 
both to equal the same entire function - thus identifying A^(z) and B_{z). The outcome is 


+ A4z)ic4z) = pA - T-{z) = 0 . 


(16) 


which requires the product factorization K{z) = Kj^{z)K-{z) and sum factorization J-{z) = 
J-+{z) B-{z) where J-{z) = F{z)/)C-{z); the product factorization is highly technical and 

is relegated to appendix the sum factorization is straightforward by inspection. After 
some algebra we obtain 

^g—i/3scos1^^ _ ^g+scosi/jj 


1 

_ COS 

except for different definitions of A^{z), B_{z) and )C+{z), IC-{z) this mirrors HK. The 
common entire function is identified, after using Liouville’s Theorem to extend to the whole 
complex plane, as zero. Thence 


_/C_(e-*^®^°®7) K,-{z)\ IC-{z) 


B-{z) ^ 

+ = 0 , 


(17) 


A+(z) 


K.+(z) 


1 


( 18 ) 


and the z-transform for the displacement coefficients B-[z) also follows. 

Technical details of the discrete Wiener-Hopf method used here are given in the appen¬ 
dices. The kernel function /C(z) is the sum of an infinite series of Hankel and modified Bessel 
functions, and its factorization is discussed in appendix together with details about the 
regularisation techniques used for evaluating various integrals. Appendix provides expla¬ 
nations and formulae required to accelerate the extremely slow convergence of the kernel’s 
series, owing to the highly oscillatory nature of the Hankel function terms and the presence of 
branch cuts. We also suggest some alternative approaches to accelerate numerical evaluation 
of the series besides the regularisation method primarily implemented here. 


3 Results 


We use two approaches to present results and illustrative examples for the displacement 
fields associated with the two types of array; in section |3.1| we consider the single semi¬ 


infinite pinned platonic grating characterised by spacing s, and in |3.2| we analyse the two- 
dimensional lattice defined by and dy. One of the approaches is the discrete Wiener-Hopf 
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method described above and we compare the results with those for a truncated semi-inhnite 
array analysed with a method attributable to [H], which we outline in section 3.1 for the 
single grating. 


3.1 The semi-infinite grating of rigid pins 


For a plane wave incident at an angle xp as in figure [^a), we determine the coefficients 
for a truncated grating by solving the algebraic system of linear equations: 

N 

'^^AkG{/3\m — k\s) = —Ui{ms), m = 1,2,N, (19) 

fc =0 

with s being the spacing of the grating’s pinned points and Uj the incident wave as defined in 
equation This is the standard Foldy scattering equation ([H]), and used amongst others 
by [ig for the Helmholtz problem, and [5] for the biharmonic plate. It is solved for N pinned 
points, and the displacements are plotted using 

N 

u{r) = Ui{r) + ^AkG{/3{\r - {ks,0)\)). (20) 

k=0 


We compare results with those obtained from the discrete Wiener-Hopf technique for the 
semi-infinite grating, which gives us the exact solution. The expressions for /C+(z), /C_(z) 


are substituted into equation (18) to determine A^(z), bearing in mind the additional singu¬ 
larity arising from the I cos i, 

term in the denominator. Referring to the expansion 
we take the inverse of the .^-transform to determine the coefficients Am- Multiplying (8) 
through by and integrating with respect to z — e*®, we obtain 


r*27r 


.2 '^A^(z) dO 



m=0 


dO = Ak2'K. 


( 21 ) 


The final result for the coefficients is the integral 

= e-^^^A+{e^^)dd, ( 22 ) 

for which the singularities and the branch cuts of appendix |B] have to be taken into account. 
The results compare well with a truncated version (at least 2000 pins) treated with the Foldy 
scattering approach. In the examples that follow we use both real parts and moduli of the 
displacement field defined by these complex coefficients, since they give us insight into the 
propagation tendencies of the scattered waves. 

For the far-field behaviour of the coefficients A^ we deduce a relation of the form 


~ 1^1 ^ 1- (23) 

The case |A| = 1 denotes the propagating Bloch wave and for |A| < 1 we observe localisation 
linked with the exponential decay of the coefficients. The ratio A is determined by 

\ — 1 ™ -iC 2 _ 

fc^oo j^A+{z)z-^dz ’ 

where the kernel in its present form is evaluated numerically. 
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3.1.1 Resonant cases 


We consider some of the frequency regimes mentioned by [5] for a finite array of pinned 
points and an infinite grating, as well as those referred to by HK for a semi-infinite grating, 
including the case they define as a resonance. This is where the incident wave travels directly 
into the end of the grating = 0 for instance) and then the diffracted wave travels parallel 
to the grating, either into (inward resonance) or away from the grating (outward resonance). 

Spectral orders of diffraction are defined according to the equation (see for instance 


0 ) 


cos 0p('0) = cos Ip + 


2ttp 

s(5 


peZ, 


(26) 


with only a finite number of the <pp{'ip) being real and representing propagating waves. The 
remaining orders are complex and represent evanescent waves. According to HK, the reso¬ 
nant cases arise when one of the spectral directions is zero (inward resonance) or vr, the case 
of outward resonance. Thus, resonances coincide with additional diffraction orders becoming 
propagating. For outward resonance 


sj3 cos 'ip -\- 27rp = —s/3, for some p = 0, ±1, ±2,..., 


(26) 


thus for p; — 0, s/3 = —wp and outward resonances occur for frequencies corresponding to /3 
being a multiple of tt. 

We observe some evidence of this effect for the case p = — 1 in figure]^ where we plot 
the scattered displacement field for values of /3 around vr for angle of incidence p) = 0 and 
period s = 1.0. We use Foldy truncated to 2000 pins and the results from Wiener-Hopf 
and Foldy are visually indistinguishable for the first several gratings, as shown in figure |^d) 
which compares the moduli of the coefficients \Ak\ in both cases for /3 = 3.1 (i.e. the black 
solid line and dots for Foldy and Wiener-Hopf respectively). 

In figure a-c) we show the real part of scattered displacement u — Ui, iov xp = 0 and 
/3 = 3.1, TT, 3.3, that is, we pass through the resonance at /I = tt. For /3 = 3.1 shown in 
figure l^a) , the high reflected energy is illustrative of the outward resonance described by 
HK for the resonant value (3 — t:. 

We contrast this with /3 — 3.3 in figure |^c) for which an additional diffraction order p = 
— 1 has become propagating. For diffraction orders p < 0 to pass off, /3s{l -\- cospx) = —27rp. 
Thus iov pj = 0, the order p = — 1 becomes propagating along with p = 0 for the resonant 
value P = TT, and the presence of two distinct propagating orders is clearly illustrated for 
P = 3.3. The two examples either side of /5 = tt show strong evidence of the circular wave 
that is typical for the end-effects of a semi-infinite scatterer. 

The resonant frequency case of /5 = tt is shown in figure |^b) and there are relatively 
lower amplitudes of displacement (with a maximum of ~ 0.3) with most of the scattering 
occurring along the grating, contrasting with the effects illustrated in figures a) and (c). 
The localisation and reduction of the maximum amplitudes to 0.3 is consistent with HK’s 
observation that for outward resonance, the cylindrical wave and all but one of the plane 
waves vanish. Indeed, the real part of the total displacement field (not shown) closely 
resembles that of the incident wave, except for the region surrounding the grating itself. 
Comparing the coefficients A^, as we do in figure |^d), further emphasises the difference 
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(a) p=3.1 (d) Coefficients 




Figure 3: Plane wave incident on a semi-infinite grating with spacing s = 1.0at'^ = 0. 
(a-c) Real part of the scattered displacement u — Ui, for /3 = 3.1, tt, 3.3. (d) Comparison of 
the moduli of the coefficients \Ak\ using Wiener-Hopf (dots for /3 = 3.1) and Foldy methods 
(shown for all). 


between this resonant case, for which \Ak\ rapidly saturates to a low value versus those for 
/3 = 3.1, 3.3 that take much longer to saturate and then do so to a larger value. 


3.1.2 Shadow boundaries 

HK refer to non-resonant cases where the diffracted field for a Helmholtz-governed semi¬ 
infinite grating consists of a set of plane waves and a cylindrical wave. The plane waves 
are consistent with those arising for the infinite grating, but do not exist everywhere. We 
obtain similar results for the biharmonic case, for which the propagating waves are due to 
the Hankel functions arising from the Helmholtz part of the Green’s function in equation 0- 
The coefficients for the scattered field are defined by equations (18) and (22) leading to the 
expression: 


1 


r*27r 


^—ikO 


dd 


27r/C_(e-*/5^‘=°W) _ ef0ei/3scosV>] ■ 

Furthermore, the representation (|^ for the displacement, leads to the following approx¬ 
imate expression for large values of 


“(C ^ XI AnHoi^Pn)- 

^ n=0 
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HK have proved in their paper that the shadow boundaries correspond to singularities 
of the ^-transform A^{z) ( |I8| ), when 2 ; is a point on the unit circle defined in the form 
2 We note that the regularised kernel ]C{z) does not have roots and poles on the 

unit circle, and hence the shadow boundaries are defined by the equation (25) and coincide 
with those discussed by HK in lEI. 


(a) p=4 



-10 0 10 20 30 40 50 


X 


(b) Coefficients 



Figure 4: The scattering of a plane wave ^ = vr/d by a grating with s = 1.0, for P = 4.0. 
(a) Real part of displacement field u. The shadow boundaries for 0o and 0_i are marked 
using dashed black lines, (b) Comparison of Wiener-Hopf (dots) and Foldy (solid black) 
coefficients \Ak\. 


Namely, denoting E{(f)) = exp{z/5scos0}, there are poles whenever F'('^) = The 

resulting values of p are labelled by and are the directions of the spectral orders of an 

infinite grating (25), of which a finite number are real and propagating. In computing the 
residues at the poles, it is merely sufficient to note that the factorization K{z) = Kj^{z)}C-{z) 
exists rather than evaluating the factors explicitly. 

The lines 9 = (ppip) act as shadow boundaries and this is illustrated in figure where 
P — 4:.0 and p — 7r/4. For p = 7r/4, the diffraction order p = —1 becomes propagating 
for P = 27r/(l + I/V 2 ) ^ 3.6806, so for P = 4.0, both p = 0 and —1 are propagating. We 
observe the two shadow boundaries for 0o(7r/4) = 7r/4 and 0_i(7r/4) = arccos (v^ — 7r)/2, 
and the field consisting of two propagating diffracted waves, as well as the reflected wave. 
Figurej^b) complements figure|^a) by showing a comparison of the moduli of the coefficients 
\Ak\ with the Foldy method represented by the solid black line, the Wiener-Hopf by the 
dots; for oblique incidence, the two methods match well and for the Wiener-Hopf we use 
the regularisation parameter = 0.0025 and 1200 intervals for the trapezoidal integration 
(see appendices . Increasing the number of intervals leads to closer convergence to the 
Foldy coefficients. 


3.1.3 Reflection and transmission 

Frequency regimes associated with total reflection and total transmission by an infinite 
grating, as identified by ISIE2], are interesting cases. For an infinite pinned grating (period 
s = 1.0, = 7r/4) the normalised reflected and transmitted energies are plotted against 

P in figure [^a). Three important values oi P = 3.2; 3.55; 3.68, which are total reflection, 
equipartition of energy and total transmission respectively, are labelled b-d and we show 
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(b) |3=3.2 


(a) 
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(C) |3=3.55 


10 20 30 

(d) |3=3.68 
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Figure 5: Plane wave incident on a semi-infinite grating with spacing s = 1.0 at an angle 
Ip = tt/A. (a) The normalised refiected, Rtot) (solid blue curve) and transmitted, Ttot, (dashed 
black) energy versus /5 for scattering by an infinite pinned grating. The real part of the 
displacement field for the semi-infinite grating is shown for (b) /5 = 3.2, (c) /5 = 3.55 and 
(d) /5 = 3.68. 


the real part of the total displacement field for the corresponding semi-infinite gratings in 
figuresgb)-(d). 

Figure [^a) shows that the maximum reflected energy (solid blue line) for the zeroth 
propagating order arises for P = 3.2. The Wood anomaly at /5 = 3.6806 signifies the passing 
off of the order p = —1, explaining the multiple orders illustrated in figure]^ for /? = 4.0. The 
real part of the total displacement field for the semi-infinite grating for the same parameter 
values p = t:/4 and /5 = 3.2 shows strong reflection to the right of the vertex, consistent 
with that observed for the full grating for the single propagating plane wave. This contrasts 
sharply with the results observed for the frequency associated with P = 4.0 in figure]^ Note 
that there is only one shadow boundary in figure [^b) (for the shadow region behind the 
grating) whereas figure |^a) shows evidence of two distinct shadow lines. 

In figure [^c) we highlight P — 3.55 which supports a mixture of reflected and trans¬ 
mitted energy for both the infinite and semi-infinite gratings; the dashed line in figure ^a) 
represents the transmitted energy. The real part of the displacement field in figure ^c) 
is consistent with the information provided by the energy diagrams for the corresponding 
infinite grating; both reflection and transmission are visible, with the shadow region above 
the grating now admitting transmittance. This transmission effect dilutes the reflection of 
figure [^b), and total reflectance is converted to total transmission by adjusting the P param¬ 
eter appropriately. This is illustrated in figure a) where P = 3.68 leads to full transmission 
for the infinite grating, and produces a similar effect for the semi-infinite grating whose total 
displacement field is shown in figure [^d) . This transmission resonance is an example of the 
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outward resonant effect for obliqne incidence. The scattered field (not shown) is reminis¬ 
cent of fignrej^b); evidence of the outward-travelling wave close to the vertex is visible in 
figure l^c) for /3 = 3.55. 


3.2 Semi-infinite lattice of rigid pins 


For an incident plane wave characterised by 'ip as in fignrej^a), we determine the coefficients 
Ak for a truncated half-plane by solving the algebraic system of linear equations (19) but 


replacing each pin with a grating in the vertical direction. We compare results with those 
obtained using discrete Wiener-Hopf. 

3.2.1 Zeros of the kernel function 


Recall from equations (12), (14) that for z = exp{iKj;dj;} lying on the nnit circle, the fnnction 
IC{z) is precisely the donbly qnasi-periodic Green’s fnnction. Therefore its zeros represent 
points on the dispersion surfaces for Bloch waves in the infinite doubly periodic structure. 
This direct connection between the kernel and the doubly periodic medium enables us to 
analyse wave phenomena nsing dispersion surfaces, band diagrams and slowness contonrs. 

We determine the scattering coefficients and displacements bk nsing the discrete 
Wiener-Hopf method. A good approximation to these results is obtained by applying Foldy’s 
method to a large enough array. For systems of at least 2000 gratings we demonstrate sev¬ 
eral interesting wave phenomena. [H] highlight the potential for varying the aspect ratio of 
rectangular lattices of pins, linking Dirac-like points with parabolic profiles in their neigh- 
bonrhood. The characteristic of this parabolic profile determines the direction of propaga¬ 
tion of localised waves. We consider first a semi-infinite rectangular lattice with aspect ratio 
7 = dy/dx = \/2. The band diagram for the doubly periodic rectangular array is shown 
in figure 7 of [E]. Here we represent the dispersion surfaces with isofrequency contour 
diagrams. 

For our initial investigations we focns on the first three dispersion snrfaces. In fignre 
we illustrate the third surface for the rectangular lattice with d^ = 1.0 and dy = \f2 with the 
Brillouin zone labelled by FXMY, where M = (7r,7r/\/2). We observe a parabolic cylinder 
profile along FX with an inflexion at F where the contours change direction for (5 ~ 5.365 
and a Dirac-like point for jd ~ 5.45 at X. We also note that the concentration of points in 
figure P tending towards the straight lines 
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Figure 6: Isofrequency contours (a selection of (5 values are labelled) for the third dispersion 
surface of the rectangular lattice with = 1, dy = \f2. The Brillouin zone FXMY is 
marked by the solid blue rectangle, with M = 


The remaining surfaces are presented together with displacement fields in the examples 
that follow. Both jd = 3.1538 for the first surface and (d = 4.40 for the second surface are 
mentioned by [19]. This immediately gives us some frequency regimes to investigate for 
finding wave phenomena including neutrality and interfacial waves. We also demonstrate 
blocking and resonant behaviour. 

3.2.2 Wave-vector diagrams 

[32] demonstrated an elegant strategy using wave-vector diagrams to investigate wave phe¬ 
nomena in planar waveguides; the technique was also outlined by [13] for photonic crystals 
in their Chapter 10. [32] was able to illustrate negative refraction, focusing and interference 
effects, having predicted them with careful analysis of the wave vector diagrams. We adopt 
an analogous approach for this platonic crystal system. 

The underlying principle is a corollary of Bloch’s theorem; in a linear system with discrete 
translational symmetry the Bloch wave vector k = (Ka,, Ky) is conserved as waves propagate, 
up to the addition of reciprocal lattice vectors. Since there is only translational symmetry 
along directions parallel to the interface, only the wave vector parallel to the interface, k||, 
is conserved. In our case, this is precisely the direction Ky. Thus for any incident plane 
wave defined by jd and k = {Kx,Ky), any reflected or refracted (transmitted) wave must 
also possess the same frequency jd'^ = uj and wave vector {K'^,i^y + 27 rl/dy) for any integer 
I and some «(.. Therefore, wave-vector diagrams may be used to analyse the reflection and 
refraction of waves within the pinned system. 

These isofrequency diagrams (also called slowness contour plots) consist of dispersion 
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curves for constant (3 (characterising the platonic crystal) and the contour for the ambient 
medium (the homogeneous biharmonic plate) on the same k^, Ky diagram (see figure]^. The 
incident wave vector (whose group velocity direction is characterised by ip) is appended to 
the ambient medium’s contour (the circle = 5.35^ = in figure [^, and a dashed 

line perpendicular to the interface of the platonic crystal and the homogeneous part of the 
biharmonic plate (k^: = 0 or FY direction here) is drawn through this point. The dashed 
line K,y = 0 and ip = 0 have been added to figure The places where this dashed line 
intersects the platonic crystal contours determine the refracted waves and their directions, 
such that their group velocity is perpendicular to the P contours, and points in the direction 
of increasing p. 

Additional information relating to the coefficients is obtained from the Wiener-Hopf 
method. Exponential decay of the form (23) with A < 1 would indicate the possibility of 
interfacial waveguide modes, with verification supported by the prediction of the refracted 
waves’ directions using wave-vector diagrams. We illustrate the technique with an intro¬ 
ductory example demonstrating a localised mode along the edge of the crystal for the first 
dispersion surface, as well as reflection and transmission which is predicted from the isofre¬ 
quency diagram. 


3.2.3 Reflection, transmission and interfacial waves for the first dispersion sur¬ 
face 

In figure [^a) we show a collection of isofrequency contours for constant P for the first 
dispersion surface for the rectangular lattice with aspect ratio = dy/d^ = \f2. This range 



Figure 7: (a) Isofrequency contours P = constant for the first dispersion surface for the 

rectangular lattice with d^ = 1.0, dy = \/2. (b) Isofrequency contours P = 3.10, 3.107, 3.11. 


of frequencies is notable for its flat bands and the parabolic profile mentioned by [H] at 
P = 3.1538 in the vicinity of a Dirac-like point. For normal incidence ip = 0, indicated by 
the direction of the arrow in figure [^a), the dashed line intersecting the ambient medium’s 
contour and normal to the interface = 0) corresponds to Ky = 0, as shown in figure [^a). 
We seek the intersections of this line with the isofrequency contour for a specific choice of 
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j3. Since the direction of the group velocity of the refracted waves is perpendicular to these 
P = constant contours, for an interfacial mode we would require the contour to be tangent 
to the line Ky = 0. 


Interfacial waves 


It is clear from figure [^a) that the contours tend towards tangency as /5 increases to /3 ~ 3.1 
but then their behaviour in the vicinity of /3 = tt is less predictable. This is connected both 


with the resonant case /3 = tt for the single semi-inhnite grating, see section 3.1.1 , and the 
parabolic profiles associated with Dirac-like points to which [11] alluded for P ~ 3.1538. This 
observation is supported by the discontinuous nature of the contours for 3.11 < P < 3.15 
when they intersect the line = 0 in figure [^a), and in particular the two distinct contour 
curves labelled by /5 = 3.15 on the right. The contours change direction at T = (0,0) for 
3.10 < P < 3.11; this corresponds to a point of inflexion on the band diagram. 

In figurej^b) we show isofrequency contours for P = 3.10 (blue), 3.107 (orange) and 3.11 
(red) to emphasise the transition of the contours through this narrow frequency window. The 
curve for /5 = 3.10 is consistent with preceding values of P but it does not touch the = 0 
line, whereas P = 3.107 (orange curve) is extremely close to touching, with the origin of a 
point of discontinuity clearly observable. This point of inflexion for some 3.107 < /5 < 3.11 
is the limit as the contours tend to = 0 and occurs at F. For this value of P, the 
refracted wave would travel along the interface, but there would be no preferential direction 
of propagation since the upper contour would indicate a “downward” direction (increasing P) 
whereas the lower contour would indicate the opposite direction. This suggests the presence 
of a standing wave, localised within the first few gratings of the half-plane of pins. 

We see evidence of this when we plot the real part of the displacement field for P = 3.107 
in figure [^b), with the predicted preferential direction of propagation indicated in figure 7(b). 
The moduli of the coefficients \Ak\ are plotted versus x-position of grating in figure 8(d), 
for the three values of P featured in figure [^b). Both P = 3.10 and P = 3.107 show 
exponential decay of the coefficients, suggesting localisation of waves, whereas the steady 
oscillation of the coefficients for P = 3.11 predicts transmission of the waves through the 
system. Figure |^b) clearly illustrates the localisation of waves within the first five gratings, 
and a clearer example of an interfacial wave is illustrated in part (a) for P = 3.10, although 
comparison with lower values of P shows similar results, but with more striking reflection 
action. 


Transmission 

Figure [^b) predicts that for P = 3.11, the preferred direction of propagation of the refracted 
waves is normal to that for P = 3.107, in the form of transmission through the system. 
Note that the opposite direction (6 = n rather than 0) is ruled out since any intersections 
corresponding to the group velocity directed towards the interface from the crystal would 
violate the boundary conditions. The coefficients’ behaviour in figure |^d) also supports 
the hypothesis of propagating waves. This is demonstrated in figure [^c) where we observe 
transmission as the wave propagates through the system without any change in direction. 
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(a) p=3.10 


(b) p=3.107 




(C) P=3.1 1 



X 


(d) Coefficients 



Figure 8: A plane wave is incident at = 0 on a lattice of 2000 gratings with = 1.0, 
dy = \f2. (a-c) Real part of total displacement field for /? = 3.10, /5 = 3.107, /? = 3.11. (d) 
Comparison of moduli of coefficients \Ak\. 


with the period of the wave’s envelope function consistent with the period for the coefficients 
in part (d). 


Coupling with finite grating stacks 

The spikes in figure |^d) for (3 = 3.10 and jd = 3.107 seem to indicate resonant interaction 
with the first few gratings of the semi-infinite array, as does the localisation evident in fig¬ 
ures |^a, b). There is a connection with the finite grating stacks analysed by [10]; resonances 
associated with the Bloch modes for the finite systems are finked to the neighbourhood of 
the Dirac-like point illustrated in figures [^a, b). 

For the pinned waveguide consisting of an odd number of gratings of period dy, a quasi- 


periodic Green’s function (12) is used to derive the dispersion equation for Bloch modes 


within the system. At each pin, the boundary conditions are u — 0. Therefore for a system 
of aligned gratings with spacing d^, 


M 


U 


(a"*) = ^ SjCj {P,md^; Ky, dy) = 0; a”* = {md^,0) 

j=-M 


for m G [—M, M], (29) 


where u is the displacement, Sj are coefficients to be determined for each Green’s function 
Gj and 2M -|- 1 is the number of gratings in the finite stack. This is equivaient to the matrix 
equation 

G S = 0, (30) 
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where S is the column vector of coefficients Sj and the Green’s function matrix G is complex 
and symmetric Toeplitz. The accompanying dispersion equation is: 

Det (G) = 0, (31) 

the solutions of which characterise the system’s Bloch modes. 


log(ldetGI) log(ldetGI) 



Figure 9: (a) Solutions of the eigenvalue problem for a system of five gratings with period 

dy = \f2 and separation d^ = 1.0. (b) Close-up of frequency window for Bloch mode 

frequencies; local maximum (3 = 3.1375 and root (3 = 3.151 marked. 


Figure [^b) features localised waves within the first five gratings so we solve the eigenvalue 
problem for a system of five gratings with period dy = \/2 and separation d^ = 1.0 in figure]^ 
The frequencies for the 5-grating system’s Bloch modes coincide with the frequency window 
for the Dirac-like point of the doubly periodic rectangular array. Similar results have been 
obtained for all 2M -|- 1-grating stacks with M < 6, with localised modes always confined 
to the range 3.10 < /I < 3.17. This is consistent with figures [^a, b) where the contours are 
well-behaved for /I < 3.10 and > 3.17, but exhibit discontinuities in the neighbourhood of 
the Dirac-like point. 

It has been shown that a finite system’s resonant Bloch modes may be coupled to appro¬ 
priately chosen incident plane waves to generate transmission action ng. We observe the 
same phenomenon for the half-plane in figure |^c) as well as evidence of similar localisation 
for l3 = 3.10, l3 = 3.107 in parts (a) and (b). Figure |^b) adds further insight when we select 
two specific values of {3 marked on the diagram; the local maximum for (3 ~ 3.1375 and the 
root (3 — 3.151 which is close to the Dirac-like point value (3 = 3.1538. In figure [l^a) the 
real parts of the coefficients Ak for /3 = 3.1538 are of order 10^, whereas for (3 — 3.1375 
in figure f^c), they are < 1. We plot the accompanying total displacement fields in fig¬ 


ure 


10 (b, d) where one should take note of the scaling bars. Both /3 values in the vicinity 


of the Dirac-like point support propagation, but the intensities of the transmitted waves are 
linked to the turning points and roots for the finite grating stack structures, illustrated in 
figure 1^ 

This coupling can be switched off by altering the period and/or separation of the grat¬ 
ings or the incoming angle of incidence. It should be noted that these examples for the 
lowest frequency band involve the zeroth propagating order only. Higher frequencies intro¬ 
duce additional diffraction orders and their coupling facilitates more exotic wave phenomena 
including negative refraction, similar to those recorded by [32] for optical waveguides. 
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Figure 10: A plane wave is incident at V' = 0 on a lattice of 2000 gratings with = 1.0, 
dy = \f2. (a) Real part of coefficients for j3 = 3.15138. (b) Real part of total displacement 
field for jd = 3.15138. (c) Real part of coefficients Ak for fd = 3.1375. (d) Real part of total 
displacement field for (d = 3.1375. 


3.2.4 Neutrality in the vicinity of Dirac-like cones 


We now consider neutrality effects that arise for parabolic profiles in the vicinity of Dirac- 
like cones. It is well known that close to Dirac points, waves may propagate as in free 
space, unaffected by any interaction with the microstructure within the crystal medium. 
These directions of neutral propagation around the Dirac point engender cloaking properties 
within the crystal, meaning that there is potential for “hiding” objects within the appropriate 
frequency regime. This property of neutrality for platonic crystals was mentioned by na, 
and was explained in a simple way in terms of the singular directions of the Green’s functions 
for the biharmonic equation. 

The first dispersion surface for the rectangular lattice with d^ = 1.0, dy = \f2 in figure]^ 
possesses such a parabolic profile near (d = tt. It is the horizontal spacing of = 1.0 that 
governs the special behaviour for normal incidence -i/i = 0 at /5 = tt, illustrated in figure 
(see equation (26)). For the vertical period dy = \/2, we observe neutral propagation of the 
incident plane wave, similar to that observed for the single grating in figure |^b), and we 
also see this for = 0 for the square lattice. Since the total displacement and incident fields 
are visually indistinguishable, we provide the real parts of the coefficients Ak in figure [T^a) 
and the real part of the scattered field in figure [T^b). It is clear that for this choice of jd = tv 
the wave does not “see” this specific lattice with d^ = 1.0, dy = v^. 
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(a) Coefficients 


X 10 ' 



(b) Scattered field for p=ji 



Figure 11 : A plane wave is incident at V' = 0 on a lattice of 4000 gratings with dx = 1 . 0 , 
dy = \f2 for /3 = TT. (a) Real parts of the coefficients A^. (b) Real part of the scattered field. 


Neutrality for oblique incidence 


For illustrative purposes we consider a square lattice of pins with aspect ratio j = dy/d^ = 1 


before moving back to the rectangular lattice with 7 = v^. In figure 12, a plane wave is 


incident at = vr/d with the Bloch parameter in the r/-direction set to Ky = /3sm{xp) = 
3.1113. 


(a) Coefficients 



(b) p=4.40, ij=:i/4 
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Figure 12 : A plane wave is incident at 'f/i = 7 r /4 on a lattice of 2000 gratings with d^ = 1.0, 
dy = 1.0 for jd — 4.40, Ky = 3.1113. (a) Modulus of the coefficients \Ak\ for both Wiener-Hopf 
(dots) with 8 = 0.005, and Foldy (solid blue). The first 50 are shown, (b) Real part of total 
displacement field. 


The real part of the total displacement field is shown in figure [I^b) , and it appears that 
the direction of the incident field is virtually unchanged. Another notable feature is that the 
stripes of the incident plane wave are replaced by spots of positive and negative intensity; 
this pattern is due to the individual scatterers and is linked to the neutrality arising from 
the Dirac-like point. Figure [T^a) illustrates the moduli of the scattering coefficients \Ak\ 
obtained using both the Wiener-Hopf (dots) and Foldy (solid blue curve) methods. The 
agreement is very good, virtually indistinguishable for the first fifty gratings, and for both 
approaches the coefficients decay to a nonzero constant supporting the neutral propagation 
that we observe in part (b). One other interesting feature is that this propagation is not the 
only wave action, the appearance of “spots” to the left of the crystal indicates that some 
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Figure 13: (a) Isofrequency contours for the first surface for square lattice with = 1.0, 

dy = 1.0. The ambient medium’s contour for /3 = 4.40 is the dashed circle + k.^. 

(b) Isofrequency curves for the 2nd surface. The dashed line Ky = 3.1113 has been drawn to 
correspond with Bloch waves in the vertical direction for (3 = 4.40, xp = it/A. The Brillouin 
zone rXMY is also shown, (c) Close-up of isofrequency contours around (3 = 4.40. (d) 
Star-like features of the dispersion curves for (3 = 4.40 at M and collection of light circles 
for the square lattice (note the multiple intersections of the light lines at Dirac-like points). 
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reflection action is also present. 

In figure[T^we use isofrequency diagrams to aid our interpretation of this neutrality effect. 
In part (a) we show the isofrequency contours for the first surface for the square lattice, and 
in part (b) the second surface since (5 = 4.40 is in the vicinity of Dirac-like cones for both. 
The contours surrounding the point M = {Kx,Ky) = (3.1113,3.1113) in the Brillouin zone 
in both parts (a) and (b) indicate that M is a Dirac-like point. The lowest value (3 contour 
near M in figure [I^a) is /? = 4.37, and in part (b) for the second surface is /? = 4.46. Part 
(c) shows a magnified picture of the neighbourhood of M, for which the (3 = 4.40 contour 
is more visible. We also illustrate an alternative derivation of the {3 = 4.40 contours in 
figure [I^d), where they take the form of star-like features at the Dirac points. Figure [I^d) 
also displays the high degree of degeneracy of the light lines at Dirac-like points such as M, 
where we observe multiple intersections. Note that figure [I^b) features the projections of 
light cone intersections onto the Ka;,Ky-plane in the form of lines ||Ka;| — = 2tt for the 

square lattice defined by 7 = 1 similar to equation ( [^ in section [3.2.1 

The dashed line Ky = = 3.1113 intersects both the ambient medium’s contour 13"^ = 

K-l + Hy and these platonic crystal isofrequency contours (3 = 4.40 precisely at the point M. 
Recalling that the group velocity should be perpendicular to these slowness contours, and 
in the direction of increasing (3, we would expect to see refracted and reflected waves at an 
angle of 7r/4 as illustrated in figure [T^b). 


3.2.5 Interfacial waves for the second dispersion surface 


but keeping all other parameter values the same. In figure [l4)(a) we plot the isofrequency 
contours for the second dispersion surface, ranging from around 3.9 to 4.65. We also add 


Now we change the aspect ratio to 7 = \/2 to obtain a rectangular lattice as in section 


3.2.3 
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(b) p=4.40 
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Figure 14: (a) Isofrequency curves for the second surface for a rectangular lattice with 

dx = dy = \/2. The angle of incidence ^ = 7r/4, the dashed line Ky = 3.1113 and 
the ambient medium’s contour = 4.4^ = ^ Ky are also shown, (b) Real part of the 

displacement field for /5 = 4.40, ijj = ii/4:^ Ky = 3.1113 for the rectangular lattice. 


the incident wave arrow for '0 = 7r/4 and the dashed line Ky = /3sin('0) = 3.1113. The 
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intersecting line at Hy = 3.1113 cuts the /3 = 4.40 contour in multiple places meaning that it 
is more dilRcult to predict the behaviour of the system, although intersections corresponding 
to the group velocity being directed towards the interface from the crystal can be ruled out 
since the only incoming energy is from the incident medium [13]. 
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Figure 15: (a) Solutions of the eigenvalue problem for a system of three gratings with 

dy = \f2 and separation = 1 for Ky = /5sin(7r/4). (b) Normalised transmitted energy 
versus jd for the triplet’s -1 diffraction order (solid black curve) and triplet’s 0 order (dashed 
red curve) for ij} = tt/A. (c) Normalised reflected energy versus jd for the triplet’s -1 diffraction 
order (solid black curve) and first pair’s 0 (dashed red curve) and -1 (thick blue) orders for 
if} = tt/A. (d) Real part of displacement for Id = 4.35, V' = 7r/4, Ky = 3.0759. 


The intersection at A indicates that the system supports transmission action in the form 
of refracted waves at an oblique angle less than 7r/4. This transmission action is visible in 
figure [T^b), although at an apparently lower intensity than the reflection observed. Another 
interesting feature of the real part of the displacement field plotted in figure [l4)(b) is a trapped 
wave close to the edge of the half-plane system. This skew-symmetric mode is localised within 
the two channels formed by the first three gratings, which is an example of wave-guiding, 
with the reflection pattern incorporating the third grating as well as the first. 

It appears that the triplet’s trapped mode is odd, so that in effect the central grating is 
not seen, and this demonstrates again the connection with finite-grating stacks discussed in 
section [3.2.3[ although here we have coupling between the diffraction orders 0 and -1. We 
explain this coupling using figure where we show the solutions to the eigenvalue problem 
for a set of three gratings with period dy = v^, separation d^ = 1.0 and Ky = /5sin(7r/4). 

Referring to figure [T^a), the first solution at /5 ~ 2.5 is due to the zeroth order only, 
whereas the solution at /3 ~ 4.35 arises for a mixture of 0 and -1 orders. This value of /3 is 
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very close to /5 = 4.40 explaining the trapped mode we observe in figure [^b), an effect we 
optimise by setting the frequency of the system with /3 = 4.35. 

We show transmitted and reflected energy profiles for the corresponding transmission 
problem for the first two and three gratings in figures 15 (b, c). In part (b), the solid black 
curve denotes the transmitted energy due to the -1 order for the triplet, explaining why the 
localised mode helps to support transmission through the rest of the half-plane system. The 
dashed red curve in part (b) represents the transmitted energy for the 0 order, emphasising 
that the -1 order dominates. This means that the first pair has to be able to facilitate 
transmission of the -1 order to reach the third grating. This is illustrated in figure [T^c) 
where reflected energy is plotted versus j3. The thick blue curve shows the reflection due to 
the pair’s -1 order is around 75% meaning that there is at least 25% transmission to feed 
the third grating, which in turn helps to feed the rest of the system. 

Figures 15 (a-c) show that the maxima for the triplets arise for /5 = 4.35 rather than 
P = 4.40. Therefore in figure [T^d) we plot the real part of the displacement field for the 
same rectangular lattice for -i/i = 7r/4, but with (5 = 4.35 and the corresponding value of Ky to 
support Bloch waves along the vertical gratings. The results are similar to those observed in 
figure [T^b) except that the transmission pattern within the pinned system is slightly more 
intense than that observed for /5 = 4.40. 


Channelling 

Recalling the second surface’s isofrequency diagram [T^a) there are several notable features 
similar to those observed for the first surface. There is a parabolic profile along FY for 
/5 ~ 4.45 and an inflexion at X for /5 ~ 4.22 where the contours change direction. Using the 
parabolic profile we are able to demonstrate an example of channelling in figure [T^ where a 
plane wave with oblique angle of incidence is bent to travel parallel to the horizontal axis. 
The wave-vector diagram is used to identify isofrequency contours parallel to = 0 for the 
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(b) p=4.40, H=:t/6 
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Figure 16: (a) Moduli (solid black) and real parts (dashed blue) of the coefficients for 

4000 gratings with /? = 4.4, -i/i = vr/G, K,y = 2.2. The first 100 of 4000 gratings with = 1, 
dy = \f2 are shown, (b) Real part of the total displacement field. 


operating frequency fd = 4.40 for the rectangular lattice with d^ = 1.0, dy = \/2; this ensures 
that the refracted waves will travel in a direction normal to the contours as illustrated by the 
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point B in figure [l^ a). The corresponding values for Ky and are 2.2 and 7r/6 respectively. 
The resultant coefficient and displacement field plots are shown in figure pT} 

Figure [l6|^a) shows very high values for both the real parts and the moduli of the 
coefficients; evidence that the system supports propagation. This is illustrated in part (b) 
where the reflection action is minimal, but the propagation through the pinned system is 
clearly visible, along with the platonic crystal’s wave-guiding effect. 


(a) Coefficients 



(b) p=4.45, n=ro'6 
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Figure 17: (a) Moduli (solid black) and real parts (dashed blue) of the coefficients Ak for 

4000 gratings with /5 = 4.45, = tt/O, Ky = 2.225. The first 30 of 4000 gratings with 

dx = \/2, dy = 1.0 are shown, (b) Real part of the total displacement field 


A similar but more intense effect is observed for (3 — 4.45 at Y with = tt/G and 
Ky = 2.225. The resultant propagation parallel to = 0 is illustrative of the uni-directional 
localised modes associated with such parabolic profiles, and can be rotated through 7r/2 by 
swapping the horizontal and vertical periods of the lattice. Thus for a lattice of 4000 gratings 
with period dy = 1.0 and spacing d^ = v^, we show an interfacial wave in figure 
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b) 

propagating in the vertical direction parallel to = 0, as predicted by the corresponding 
isofrequency contour, and very little wave action inside the pins as predicted by the decay 
of the coefficients in figure [l7|(a) . 


3.2.6 Higher frequencies - interfacial waves 

It is clear that as the frequency oj = (3“^ increases, a more complicated picture emerges. 
The third surface illustrated in figure and magnified here in figure [I^c) possesses many 
important features. At the point F the contours change direction by 7r/2 for /5 ~ 5.365 
and there is a Dirac-like point at A for /5 ~ 5.45. This narrow frequency window supports 
an array of interesting wave phenomena, similar to those discussed for the first surface in 
section 13.2.31 

We begin by considering a neighbouring frequency (3 = 5.6 illustrated in figure [T^c) 
where we show a collection of isofrequency contours for the platonic crystal’s third dispersion 
surface, as well as the contour for the homogeneous biharmonic plate medium - an arc of 
the dashed circle defined by (3'^ 


K^ + i^y for the case [3 = 5.6. We also refer again to 
equation (28) to explain the appearance of the light cone projections on the isofrequency 
diagram. 
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(a) Coefficients (b) 6=5.6, -u=0.1361 



Figure 18: A plane wave is incident at ^ = 0.1361 on a lattice of 2000 gratings with 
dx = 1.0, dy = \/2 for /? = 5.6, Ky = 0.76. (a) Moduli of the coefficients \Ak\ with the first 
50 shown for both Foldy (solid blue curve) and Wiener-Hopf (dots), (b) Real part of the 
total displacement field, (c) Magnification of the third dispersion surface, (d) Real part of 
the total displacement field for jd = 5.6, = 0.2712, Ky = 1.5. 


We seek an interfacial wave so we choose ip such that Ky = [3 sin('^) is tangent to the 
/5 = 5.6 contour for the third surface (point A in figure [I^c)). The vertical arrow predicts the 
propagation direction of the refracted waves for this Ky = 0.76 and corresponding ip — 0.1361. 
In this way we predict that some refracted waves will propagate perpendicularly to this 
contour, i.e. in a vertical direction parallel to the interface. We plot the real part of the 
displacement field for these parameter settings in figure [l8|(b) , along with the moduli for the 
coefficients \Ak\ in figure [T^a) generated by both the Foldy and Wiener-Hopf methods. 

Figure [l8|(b) shows low propagation within the pinned structure, and relatively low re¬ 
flection, compared with figure [T^b) for example. However there is clearly a travelling wave 
localised within the first few gratings, albeit with a relatively long wavelength. The expo¬ 
nential decay of the moduli of the scattering coefficients \Ak\ is illustrated in figure [T^a). 
In figure [I^d) we consider an example for an arbitrary value of Ky to highlight the contrast 
in behaviour for various ip and Ky. For Ky = 1.5 and ip = 0.2712, we show the real part 
of the total displacement field in figure 18 (d). Although there is increased intensity at the 
interface, as one would expect, the dominant behaviour is propagation and in almost the 
same direction as the incident wave. This is predicted by the direction of the blue arrows at 
B and on the /3 = 5.6 contour close to the incident wave’s arrow ip = 0.2712 in figure (c). 

Finally we consider the case ip = 0, Ky = 0. At the point F in figure [T^c), the contours 
for (3 = 5.35 support propagation in a direction normal to those of /3 = 5.40, which support 
propagation parallel to the interface similar to A for /3 = 5.6 for ip = 0.1361 in figures [I^b,c). 


26 




































(a) p=5.40 


(b) p=5.45 
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(c) Coefficients 



Figure 19: A plane wave is incident at V' = 0 on a lattice of 4000 gratings with dx — 1.0, 
dy = \f2. (a, b) Real part of total displacement field for /3 = 5.40, /3 = 5.45. (c) Comparison 
of the moduli of coefficients \Ak\. 


Thus there is a point of inflexion for 5.35 < (3 < 5.40 and it occurs for (3 ~ 5.365. For 
(3 < 5.365, the refracted waves propagate through the pinned system rather than along its 
edge, which is what happens for 5.365 < (3 < 5.45. However at the Dirac-like point X; 
P ~ 5.45, the waves exhibit a mixture of edge localisation and propagation through the 
system. This is because we have the third, fourth and fifth dispersion surfaces meeting at 
this frequency, leading to a combination of behaviours. The third surface [I^c) predicts 
propagation along the edge for /3 = 5.45 contours while the fourth and fifth surfaces predict 
propagation parallel to K,y = 0, into the pinned lattice. 

We observe these properties in figure [I^ where parts (a,b) show the real part of the total 
displacement field for P = 5.40 and P = 5.45. The amplitudes are extraordinarily large 
for P = 5.365 where the isofrequency contours are about to change direction to support 
interfacial waves rather than propagation through the pinned system (see F in figure [^c)). 
This transition has occurred in figure [T^ a) for P = 5.40 where the localised interface mode 
clearly dominates any action inside the pinned region. The comparison of the coefficients 
Ak in figure [I^c) also demonstrates the decay of the coefficients to zero inside the crystal 
for P = 5.40. However P = 5.45 close to a Dirac-like point supports both edge localisation 
and wave propagation as illustrated in figure 19 (b, c), where a wave is seen to propagate 
through the pinned lattice with the wavelength of the envelope function matching that of 
the coefficients. 
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4 Concluding remarks - interfacial waves and dynamic 
localisation 

The theory and examples presented in this paper have identified novel regimes, typical of 
flexural waves in structured Kirchhoff-Love plates, for the case of a semi-infinite structured 
array of rigid pins in an otherwise homogeneous plate. The interface waveguide modes have 
been identified and studied here. 

For a half-plane occupied by periodically distributed rigid pins, we have demonstrated 
an interplay between the transmission/reflection properties at the interface and dispersion 
properties of Floquet-Bloch waves in an infinite doubly periodic constrained plate. Specifi¬ 
cally, we have analysed regimes corresponding to frequencies and wave vectors that determine 
stationary points on the dispersion surfaces, as well as Dirac-like points. For such regimes 
we have demonstrated that the structure supports localised interfacial waves, amongst other 
dynamic effects. 

The localisation was predicted by an analytical solution, and a formal connection has 
been shown here between the doubly quasi-periodic Green’s function for an infinite plane 
and the system of equations required to obtain the intensities of sources at the rigid pins, 
which occupy the half-plane. 
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A Factorizing lC{z) 


A critical technical detail is the factorization K,{z) = K,j^{z)K,-{z) with /C+,/C_ respectively 
analytic inside and outside the unit circle. Laurent’s theorem gives 



(Al) 


where G+ is a circle of radius c_|_ slightly larger than the unit circle, and G_ is a circle of 
radius c_ slightly smaller than the unit circle (see figure ; the desired factorization is 



K+{z)K-{z). (A2) 
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We note here that g is of the form for (7+ and C-, with z = 

that C± are chosen so that log K, is well defined. 


Referring to the factorization (A2), 




log/C(£)) 


c± Q-z 


dg 


It is also assumed 


is evaluated for z = for some argument 9i on the unit circle, with g = c±e*^ for 0 < 0 < 
27r. Here we outline two possible choices of regularisation, both of which give good results. 
Defining g = with g± = exp{±h}, we have 

dg = ig±e^^d6 


and 


Then we may write 


JC±{z) = exp 


±1 log/C(p±e*^)p±e*^ 

2tt Jo g±e^ - z 


JC±{z) = exp 


^ logJC{g±e^^) 1 
27r Jo 1 ~ j 


exp 


^ r ■°g^<pf2 d4. 

27r Jo Q± - / 


(A3) 


The alternative way is to define /ds = P + iS and then determine c+,c_ accordingly. Both 
approaches achieve good and very similar results (comparable to order 10“®). 


B Accelerated convergence for JC{z) 

The biharmonic operator’s kernel /C(z) converges extremely slowly because of the highly 
oscillatory nature of the Hankel function terms. The direct correspondence with the quasi- 
periodic grating Green’s function for a specific definition of z = exp{ma;s} enables us to 
implement some known accelerated convergence techniques. Twersky [28] investigated the 
convergence of the Schlomilch series 

OO 

Z 21 ipD) cos{pD sin po ), (Bl) 

p=l 

OO 

Z 21+1 (pD) sm{pD sin po) , (B2) 

p=l 

where is the nth order Bessel function and D > 0, 0 < sin^o < 1- In particular, we 
consider the series in the form 

OO 

Tin = '^H^\pD)[exp{ipDsmi/jo}{-l)'^ + exp{-ip D sin po}]. (B3) 

p=i 

Here Hn^ = + iA„ is the standard Hankel function of the first kind, with Jn and A„ being 

the Bessel and Neumann functions. It is clear that D can be replaced by /5s to coincide with 
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our treatment, and we note that /5s sin(';/^o) is associated with the Bloch parameter for 
the periodic structure we are considering. 


It is well-known that the representation (|B3| is too slowly convergent for practical use. 
Twer sky 


derived an alternative rapidly convergent representation in terms of elementary 
functions. However he did not evaluate the n = 0 case which we require here, but instead 
referred to the work of ng. Rewriting the representation for the kernel function (2.9) in the 
form 

■ 9' 

= + -:^Ko{/3sj)] {z^ + ^"9 + (B4) 


we substitute n = 0 in (B3) to obtain the series 


2 (pD) cos {p D sin xpo), 

p=i 


(B6) 


which is immediately associated with the slowly convergent part of representation (B4) for 
2 ; on the unit circle, with D = f3s and 2 ; = 6 = /5s sin= i^xS- For 2 ; sitting precisely 

on the unit circle, we may use the accelerated convergence formulae of [22] which involve 


(B6) 


grating sums Sq for the Hankel functions, and Sq for the Bessel K functions: 


s « + +1 = ? 

TT S 


E 


1 

Xp 


1 

Xp 


'E 


1 


\Xp 


1 

Xp J . 


where the right-hand sum is made up of propagating and evanescent parts, and is cubically 
convergent. The terms Xp Xp = iTp, Tp > 0 are dehned by 


Kjp 

27rp 
l^x + ) 

s 


(B7) 

Xp 

/ y/j" - 

1 VV “ 

, Kp^ > /5^ 

(B8) 

= 

7/5^ + 


(B9) 


where p G Z. There is a finite number of propagating orders p (the 1/Xp terms in (B6)) with 
all other orders being evanescent. 

Branch cuts arise from the Helmholtz part of the kernel function K(z) (2.9) and not 
from the modified Helmholtz part. This is illustrated in figure |B1[ where part (a) shows 
the contribution from the sum of Hankel functions, and part (b) shows the real part of the 
Bessel K contributions (the imaginary part is of order 10“^°). Not only does figure |BT| (b) 
illustrate the exponentially small contribution from the modified Helmholtz operator, but 
also that it is a well-behaved function with no branch cuts. These branch cuts arise from 
the grating sums which include a factor of the form l/(/5^ — + 27rp/s)^), which is often 

used to define the “light lines” for the system. 
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(a) 


(b) 




Figure Bl: (a) Imaginary part of the Helmholtz part of the kernel function for P = 4.0, 
s = 1.0 with = e*® for 0 G [3.95,4.05] using 5000 terms (solid grey line) and convergent 
grating sums (dashed blue line), (b) Real part of the modified Helmholtz part of the kernel 
function for the same parameter values with 5000 terms. 


B.l Series representations for K,{z) for numerical evaluation 

The accelerated convergence formulae of equation ( |B6[ ) are only valid for lC{z) if 2 ; lies 
on the unit circle. Expressions for }C+{z) and K.-{z) in (A2) involve g which lies either 
just inside or outside the unit circle. Therefore to evaluate /C+ and /C_ we must use the 
direct form for )C{z) (|B4[) with either regularisation or the use of a remainder function which 


employs a finite number of terms directly and adds an infinite tail evaluated using asymptotic 
approximations. For the latter treatment, we rewrite IC(z) in the form 


E 


^ j=N+l 


H^o\M) + -Ko{psj) 

TT 


2i 


N 


+ — T 

8/52 ^ 


Ho\Psj) + —KoiPsj) 

TT 


{z^ + z 




(BIO) 


where N denotes a finite number of terms for direct application of the kernel series, with 
the remainder evaluated via a function we define as R{z) based on the asymptotic analysis 
by HK. We write 


00 

R{z)= XI + ^~^) 

i=iv+l 


n=l 


(Bll) 


where we have used the change of index of summation n = j — N. 

We replace the Hankel functions by their asymptotic forms for large Ps, since we are only 
considering the tail oi IC{z). Hence, 


R{z) = e 


./4 ^ {ze^^^Y+^ 

V /5s ^ V7r(n +W) 


+ e' 


74 /Zv_(if£l 

V /5s ^ A/7r(n +A)’ 


(B12) 
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for which we define the function F[z) by 


m = E 


yU + N 


^ y7r(n + N) tt ^ Jo 


OO ^oo 

V / (B13) 

_ 1 Jo 


where we refer to Appendix 1 of HK. For \z\ < 1, we may interchange the order of summation 
and integration such that 


F{z} = 


^ «Qo 

- b = ze-*\ | 6 | < 1 . 

^Jo ti 


Thus, 


nz) = l[ 

^ Jo 


2 6^+1 dt 


i-b 


and 
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F{ze^^^) + F{ 
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Thus, referring to equation (BIO) we obtain 


(B14) 


(B15) 


IC{z) = 


Sffi 


N r 


E 

j = l L 


+ -Ko{/3sj) 


2i 


TT 


{z^ F z ^) + R{z) + 1 


(B16) 


This more amenable representation for K{z) enables us to determine explicit expressions 
for /C+(z) and JC-{z) of equation (A2) which are required to evaluate A+(z) (2.18). However 
since regularisation is our preferred method of handling the branch cuts, the introduction 
oi 15s = /J + iS ensures that the sums converge sufficiently quickly to avoid the necessity of 
using the remainder function R{z) for }C+{z) and K,-{z) away from the unit circle. 
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